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We examine the performance of various commonly used integration schemes in dissipative particle dynamics 
simulations. We consider this issue using three different model systems, which characterize a variety of differ- 
ent conditions often studied in simulations. Specifically we clarify the performance of integration schemes in 
hybrid models, which combine microscopic and meso-scale descriptions of different particles using both soft 
and hard interactions. We find that in all three model systems many commonly used integrators may give rise to 
surprisingly pronounced artifacts in physical observables such as the radial distribution function, the compress- 
ibility, and the tracer diffusion coefficient. The artifacts are found to be strongest in systems, where interparticle 
interactions are soft and predominated by random and dissipative forces, while in systems governed by con- 
servative interactions the artifacts are weaker. Our results suggest that the quality of any integration scheme 
employed is crucial in all cases where the role of random and dissipative forces is important, including hybrid 
models where the solvent is described in terms of soft potentials. Regarding the integration schemes, the best 
overall performance is found for integrators in which the velocity dependence of dissipative forces is taken into 
account, and particularly good performance is found for an approach in which velocities and dissipative forces 
are determined self-consistently. Remaining temperature deviations from the desired limit can be corrected by 
carrying out the self-consistent integration in conjunction with an auxiliary thermostat, in a manner that is simi- 
lar in spirit to the well-known Nose-Hoover thermostat. Further, we show that conservative interactions can play 
a significant role in describing the transport properties of simple fluids, in contrast to approximations often made 
in deriving analytical theories. In general, our results illustrate the main problems associated with simulation 
methods in which dissipative forces are velocity dependent, and point to the need to develop new techniques to 
resolve these issues. 



I. INTRODUCTION 



One of the greatest challenges in theoretical physics is to 
understand the basic principles that govern soft condensed 
matter systems, such as polymer solutions and melts, colloidal 
suspensions, and various biological processes. Experimental 
studies of these complex systems are often complemented by 
numerical simulations of model systems, which can provide 
a great deal of information not easily accessible by experi- 
ment. In this regard, molecular dynamics |jl]] (MD) is often the 
method of choice, and indeed it can elucidate various physical 
phenomena on a microscopic level. In general, however, such 
an atomistic approach is problematic since many intriguing 
processes in soft matter systems are not dictated by micro- 
scopic details but rather take place at mesoscopic length and 
time scales (roughly 1-1000 nm and 1-1000 ns) which are be- 
yond the practical limits of MD. In such cases, it is necessary 
to model soft matter systems by viewing them from a larger 
perspective than from a microscopic point of view. In practical 



terms, this means that one has to design ways to simplify the 
underlying systems as much as possible, while still retaining 
the key properties which are expected to govern the processes 
of interest. Recently, this approach has attracted wider atten- 
tion as various "coarse-grained" simulation techniques have 
been developed [| f| |, [| ^ ||] with the purpose of study- 
ing mesoscopic physical properties of model systems. 

Dissipative particle dynamics ^, [7|, ||] (DPD) is partic- 
ularly well suited for this purpose. DPD is characterized by 
coarse-graining in particle representation, which allows stud- 
ies of systems at mesoscopic length scales and a simplified 
description of inter-particle interactions allowing for stud- 
ies at mesoscopic time scales. Since DPD preserves hydrody- 
namic modes, one may characterize DPD as momentum con- 
serving Brownian dynamics. For these reasons, DPD is a very 
promising method for mesoscopic studies of soft systems and 
recently has attracted considerable interest in studies of poly- 
mers, microphase separation, ^ and lipid bilayers, jl0| ] 
among others. 

Despite its advantages, DPD has certain practical problems 
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that have to be resolved before extensive use in large-scale 
simulations. Many of them are related to the idea of coarse 
graining which can be done by simplifying molecular repre- 
sentations, and then replacing the "fast" variables related to 
the coarse-grained degrees of freedom by random noise. The 
random noise mimics thermal fluctuations and hence drives 
the system. In DPD, this idea is implemented by a special 
"DPD thermostat" 0, ||| in terms of dissipative as well as 
random pairwise forces such that the momentum is conserved. 
This is a prerequisite for the emergence of hydrodynamic flow 
effects on a macroscopic scale. However, due to the DPD ther- 
mostat and the resulting stochastic nature of the equations of 
motion, the quest for a suitable integration scheme in DPD is a 
non-trivial task. It has been recently observed that various in- 
tegration schemes commonly used in classical MD lead to dis- 
tinct deviations from the true equilibrium behavior, including 
deviations from the temperature predicted by the fluctuation- 
dissipation theorem, and artificial structures observed in the 
radial distribution function. [0, |IJ [[J, [jj, These findings 
demonstrate the serious practical problems associated with the 
use of DPD and raise concerns regarding its future application 
to large-scale simulations of soft systems. 

A related problem regards hybrid models, where the aim is 
to combine microscopic models of biomolecules with a meso- 
scale modeling of the solvent. [15 [Tq] In this promising 



approach, one can examine microscopic properties of com- 
plex biological molecules in an explicit solvent but with a re- 
duced computational cost. While biomolecules are described 
by hard conservative interactions such as Lennard-Jones and 
Coulombic forces, the solvent can be described by DPD as 
a softly interacting fluid. The drawback is that the integra- 
tion schemes may again lead to deviations from the true equi- 
librium behavior. To our knowledge, the role of integration 
schemes in these cases, where both soft and hard interactions 
are used within a meso-scale DPD simulation, has not been 
studied yet. These examples clearly highlight the current need 
to examine the relative performance of different integration 
schemes in DPD under various conditions, and develop new 
integration techniques where the special features of DPD are 
properly accounted for. 

In this work, we study the performance of a number of com- 
monly used integration schemes in DPD simulations. They 
all are based on the velocity- Verlet scheme but differ in how 
the velocity dependence of dissipative forces in DPD is taken 
into account. We test the integrators by studying a number of 
physical observables such as the temperature, the compress- 
ibility and the tracer diffusion coefficient, and evaluate their 
performance in three different model systems. We first ex- 
amine how the integrators perform in the absence of conser- 
vative forces. This case was partly discussed in our previ- 
ous work, $1% which is here extended by a thorough discus- 
sion of the results and the details of self-consistent integration 
schemes suggested in Ref. [n|. Then, by increasing the rela- 
tive importance of conservative forces, we eventually obtain a 
model which is used to assess the performance of integration 
schemes in a hybrid approach. 

We find that various commonly used integration schemes 
in MD and DPD indeed lead to pronounced artifacts in actual 



physical quantities. These artifacts are found to be strongest 
in weakly interacting systems, where interactions are soft and 
dominated by random and dissipative forces. In the opposite 
limit, where hard conservative interactions govern the system 
under study, the artifacts due to integration schemes are less 
pronounced. We conclude that the quality of an integration 
scheme employed is crucial in all cases where the role of 
random and dissipative forces is important, including hybrid 
models where the solvent is described in terms of soft poten- 
tials. 

Regarding the integration schemes, best overall perfor- 
mance is found for an integration scheme which involves the 
self-consistent determination of particle velocities and dissi- 
pative forces. For cases where precise temperature control is 
crucial, we further suggest and analyze in detail an additional 
auxiliary thermostat which corrects for the residual tempera- 
ture deviations. 

The outline of the paper is as follows. In Sect. 0, we first 
review the essential background of DPD and then introduce 
the three model systems studied in this work. The integra- 
tors which are tested are described in Sect, |5|, after which, in 
Sect. |[v|, we present and discuss the test results. In Sect. [v| we 
discuss the special case of tracer diffusion behavior of DPD 
particles and compare our findings to previous theoretical de- 
scriptions. Finally, we close this paper with a short summary 
and discussion in Sect. fVll 



II. METHODS AND MODELS 

Below we give a short summary of DPD and describe the 
model systems used in this study. For more thorough accounts 
on DPD, see Refs. | and |. 



A. Dissipative Particle Dynamics 

In the present work, we study a simple model fluid sys- 
tem described by N identical particles each with mass m, and 
which have coordinates ri, and velocities ili. Interparticle in- 
teractions are characterized by the pairwise conservative, dis- 
sipative, and random forces exerted on particle "i" by particle 
"j", respectively, and are given by 
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where r 



and v 



j ' 



The are symmetric random variables with zero 
mean and unit variance, and are independent for different pairs 
of particles and different times. 

The pairwise conservative force of Eq. (jl]) is written in 
terms of a weight function w c {r^ ■), whose choice is dictated 
by the system under study. In principle, u) G {rij) can include 
various kinds of forces due to e.g. electrostatic interactions, 
as well as descriptions of detailed intermolecular interactions 
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such as van der Waals forces. However, since DPD has been 
designed to model molecular systems on a mesoscopic level, 
a detailed atomistic description of interactions is, in many 
cases, not necessary. Instead, it is often preferable to use 
soft-repulsive interactions of a relatively simple form. This 
approach is justified by observations by Forrest and Suter that 
coarse graining of a molecular representation tends to soften 
interactions. [ |l8| | Recent work by Flekk0y et al. [ |l9| ] also 



supports this view. We will return to this issue in Sect. [I B 
where the actual form of lo c (fy) will be discussed in more 
detail. 

Unlike the conservative force, the weight functions lo d (r^ ) 
and w R (rij) of the dissipative and random forces cannot be 
chosen independently. Physically, and F3 have to be cou- 
pled, since thermal heat generated by the random force must 
be balanced by dissipation. The precise relationship between 
these two forces is determined by the fluctuation-dissipative 
theorem, which sets conditions for both the weight functions 

w D (r ij ) = [u a (r iJ )] 2 (4) 

and the amplitudes of these forces 

a 2 = 2-/k B T* , (5) 

where T* is the canonical temperature of the system. 

DPD samples phase space according to the canonical en- 
semble with a potential energy determined by the conserva- 
tive force described by Eq. djl). Consequently, static proper- 
ties such as the pair correlation function and the specific heat 
could be calculated equally well using any stochastic tech- 
nique (such as the off-lattice Monte Carlo method). On the 
other hand, as in molecular dynamics, DPD also provides a 
means to calculate dynamical properties of the system. Con- 
sequently, a method is required to evolve the system in time, 
which in DPD is usually done by integrating Newton's equa- 
tions of motion. Unlike the case in standard molecular dy- 
namics, the presence of a stochastic contribution to the force 
in DPD implies that the equations of motion are now given by 
the set of stochastic differential equations 



dn = 
dVi = 



Hi dt , 

J-fPfdt + F^dt 
m, \ 



FrVdt 



(6) 
(7) 



where F R = V . Fff is the total random force acting on 

particle "£", and Ff and Ff are defined correspondingly. The 
velocity increment due to the random force in Eq. (0) is writ- 
ten in a form which can be given a precise meaning by iden- 
tifying it as the infinitesimal increment of a Wiener process. 
PQ ] In practice, finite time increments are used in the simula- 
tions, and the equations of motion [Eqs. (||) and (Q)] have to 
be solved by some integration procedure. As will be seen in 
Sect. £v|, this may lead to serious artifacts if the special fea- 
tures of DPD are not taken into account. 



B. Model systems 

In this study, we investigate the performance of various in- 
tegration schemes using three different model systems. They 



all are based on a 3D simple model fluid system with a fixed 
number of similar spherical particles. The differences be- 
tween the model systems arise from interaction effects, which 
are varied step by step from an ideal gas to a more realistic 
description of an interacting fluid system. The models studied 
here are described below. 



1. Model A 

We first consider the case characterized by the absence of 
conservative forces (a = 0). This choice corresponds to an 
ideal gas (sometimes termed "ideal DPD fluid" within the 
framework of DPD), which provides us with some exact the- 
oretical results to be compared with those of model simula- 
tions. Here, the dynamics of the system arise only from ran- 
dom noise and the dissipative coupling between pairs of parti- 
cles. The random force strength is chosen as a = 3 in units of 
ksT*, and the strength of the dissipative force 7 is then deter- 
mined by the fluctuation-dissipation relation in Eq. (|j). The 
weight function uo R (rij) was chosen as in various previous 
works, [0 



uj R {n 3 ) 



1 - r tj /r c , for r tj < r c 
, for r tj > r c 



(8) 



where r c is a cut-off distance. The weight function ui D {rij) 
is defined via Eq. (|j). Therefore, the dissipative and random 
forces are just soft pairwise repulsions acting along the line of 
centers of two DPD particles, and 7 and a are the amplitudes 
which define the maxima of these forces. 

In our simulations we use a 3D box of size 10 x 10 x 10 
with periodic boundary conditions, where the length scale is 
defined by r c = 1, and a particle number density p = 4. 

2. Model B 

Model B is a simple interacting DPD fluid. Its main differ- 
ence with respect to model A is the presence of a conserva- 
tive force, which we choose to have a strength a = 25 and a 
weight function u) 1 (rjj) of the same form as the random force 
in Eq. (|J). In all other respects, this model system is identical 
to model A. 



3. Model C 

Model C is a variation of model B. Instead of soft poten- 
tials, we now use hard conservative interactions. The conser- 
vative potential between particles "i" and "j" is described by 
the truncated and shifted Lennard- Jones potential 




1 ^ 



(9) 

such that the potential is purely repulsive and decays smoothly 
to zero at r c . We choose I — and e = ksT*, and 
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therefore r r 



£2 1 / 6 = 1. Forr^ > r c , C/g(r tf ) 



0. 

?c _ 



The pairwise conservative force follows directly from = 
—^7UH. The dissipative and random forces are described by 
Eq. $>. 

Simulations were carried out in a 3D box of size 16 x 16 x 16 
with periodic boundary conditions, with a ranging from 1 to 
200, and with particle densities p = 0.1 and p = 0.7. 

Finally, let us briefly justify the choice of these model sys- 
tems. In model A, the idea is to study integrator-induced 
artifacts in a case, where the role of random noise with re- 
spect to conservative interactions is as pronounced as possible. 
Model B corresponds to a typical situation where large-scale 
processes such as microphase separation and morphological 
properties of complex systems are studied in terms of coarse- 
grained particles. In that case, details of molecular represen- 
tation are no longer accounted for, and all interactions are de- 
scribed in terms of soft potentials. Finally, model C aims to 
gauge integrator-induced effects in a hybrid approach, where 
both hard and soft interactions are present. One likely scenario 
of this idea is to model solute molecules with realistic atomic 
force fields, while the solvent is described in a coarse grained 
fashion. In the present work we restrict our test simulations 
to simple spherical particles which interact via Lennard-Jones 
type interactions, since that should already allow a reliable as- 
sessment of integrator-induced artifacts in hybrid models. 



C. Practical details 

One of most important practical aspects within DPD is the 
stochastic nature of the interactions. This is built in to the ran- 
dom force of Eq. (Q) via £y, which are independent random 
variables with zero mean and unit variance. In the present 
work, we have described them by uniformly distributed ran- 
dom numbers u € U(0,1) such that, for every pair of par- 
ticles at any moment, we generate a different stochastic term 
£ = \/3(2u — 1). This approach is very efficient and yields re- 
sults that are indistinguishable from those generated by Gaus- 
sian random numbers. [ pT| 

In generating the random numbers, we used a pseudoran- 
dom number generator RAN2, which is based on the 32-bit 
combination generator first proposed by L'Ecuyer [22| and 
later published in Numerical Recipes [ [23[ | using shuffling. In 
a recent study, p4] ] where several pseudorandom number gen- 
erators were tested in DPD model simulations, it was demon- 
strated that RAN2 performs very well in simulations of simple 
fluids. 

The length scale in the simulations is defined by r c = 1 and 



the time scale is given in units of r c y/mJkgT* . The energy 
scale is defined by setting the desired thermal energy to unity 
via k B T* = 1. 

The simulations were carried out with particle numbers of 
the order of a few thousand (4000 in models A and B, and 
roughly 2800 in model C for p = 0.7). The number of time 
steps varied depending on the size of the time increment At 
such that the total simulation time was about 5000 - 10000 (in 
units of r c ^/m/fcsT*). 



III. INTEGRATORS 

One of the central issues in molecular dynamics calcula- 
tions is the integration of the equations of motion. In the con- 
text of MD, the present understanding of this issue is rather 
clear and comprehensive. J25| ] However, in the case of DPD 
simulations the situation is more problematic. To clarify this, 
let us consider the equations of motion in detail. Using Eq. (0) 
for the velocity term we obtain 



dvi = — 

mi 



■dt \rij) e K 



dt ^lo 1 

3& 



which immediately reveals two potential problems. First, due 
to the stochastic nature of interactions, the time reversibil- 
ity is no longer guaranteed. Another serious problem arises 
from the dissipative forces, which depend on the pairwise ve- 
locities of all pairs of particles. This seemingly minor de- 
tail is absent from classical MD simulations but leads to sig- 
nificant problems in DPD simulations, including artifacts in 
various physical quantities measured from simulation studies. 
|Q, [jr], [l2| |3]| In principle, this problem could be solved by 
finding a self-consistent solution for both dissipative forces 
and particle velocities by inverting an appropriate interaction 
matrix of size A x A at every time step. However, it is obvi- 
ous that this approach is generally not feasible, and thus one 
must search for more practical solutions. 



A. Simple velocity- Verlet based integration schemes 



We use the velocity-Verlet scheme |J2(J as a starting point 
and consider various previously used integrators based on 
this approach. These are summarized in Table [| where the 
acronym "MD-VV" corresponds to the standard velocity- 
Verlet algorithm used in classical MD simulations. The MD- 
VV scheme is (in the case of solely conservative forces) 
a time-reversible and symplectic second-order integration 
scheme, which has been shown to be relatively accurate in 
typical MD simulations especially at large time steps. [^7]| Al- 
though some higher-order algorithms are more accurate than 
the MD-VV, iQ the simplicity of MD-VV makes it a good 
starting point for further development. 

Unlike in molecular dynamics, the forces in DPD depend 
on the velocities. This fact is not accounted for within the 
MD-VV scheme. In an attempt to deal with this complica- 
tion, Groot and Warren subsequently proposed ^ a modified 
velocity-Verlet integrator ["GW(A)" in Table j]. In this ap- 
proach, the forces are still updated only once per integration 
step, but the dissipative forces are evaluated based on inter- 
mediate "predicted" velocities v.° . The underlying idea of v° 
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is to use a phenomenological tuning parameter A, which mim- 
ics higher-order corrections in the integration procedure. The 
case A = 1/2 corresponds to the usual MD-V V, while other 
suggested choices range from zero to one. The problem is 
that the relative merits of different numerical values for A are 
poorly understood. Groot and Warren studied the tempera- 
ture control in an interacting fluid and found that A = 0.65 
works better than A = 1/2. [Q| In a different study, Novik 
and Coveney concluded that A = 1/2 gives a more accurate 
temperature than A = 1. [ pi] ] Thus, it is evident that the op- 
timal value of A, which minimizes the temperature shift and 
other artifacts, depends on model parameters and has to be 
determined empirically. 

Recently, Gibson et al. proposed Jl3[ | a slightly modified 
version of the GW integrator. This "GCC(A)" integrator up- 
dates the dissipative forces [step (5) in Table |[ for a second 
time at the end of each integration step. This approach suffers 
from the same problem as the GW integration scheme, i.e., it 
uses a phenomenological parameter whose optimization de- 
pends on the system and the conditions that are being mod- 
eled. Based on a few model studies by Gibson et al., values 
of A between 1/2 and 1 may be preferable to smaller values. 



Despite the use of a phenomenological parameter, the GCC 
scheme is a promising approach for DPD simulations. A ratio- 
nal starting point is to fix A to a value of 1/2, which leads to an 
integrator equivalent to the MD-V V scheme supplemented by 
the second update of the dissipative forces. This Verlet-type 
integrator, here termed "DPD-VV", is particularly appealing 
because it does not involve any tuning parameters, yet it takes 
the velocity-dependence of the dissipative forces at least ap- 
proximately into account. In addition, it is computationally 
very efficient since the additional update of dissipative forces 
is an easy task compared to the time-consuming part of updat- 
ing neighbor tables. 

In this work, we consider besides the schemes GW(A = 
1/2) = MD-VV and GCC(A = 1/2) = DPD-VV also 
GW(A = 0.65) and GCC(A = 0.65). 



B. Self-consistent velocity- Verlet integrators 



Unfortunately, as will be shown in Sect. IV, all of the above 
integrators display pronounced unphysical artifacts in the ra- 
dial distribution function g(r), and thus do not produce the 
correct equilibrium properties (see results and discussion be- 
low). This highlights the need for an approach in which the 
velocity-dependence of dissipative forces is fully taken into 
account. In principle this problem can be easily addressed 
by solving the velocities and dissipative forces in a self- 
consistent fashion. In practice, however, there is no unique 
way to do this. In this work, we present in Table || the update 
schemes for two self-consistent schemes which are variants of 
DPD-VV. The basic variant SC-VV, which is similar in spirit 
to the self-consistent leap-frog scheme introduced by Pago- 
nabarraga et al, [|l2|| determines the velocities and dissipative 
forces self-consistently through functional iteration, and the 
convergence of the iteration process is monitored by the in- 



stantaneous temperature fc^T. 

In the second approach, which we call SC-Th, we couple 
the system to an auxiliary thermostat and obtain an "extended- 
system" method in the spirit of Nose-Hoover. [ 29 1 The idea 
behind this approach is that whenever (ksT) deviates from 
ksT*, the dissipation rate is on average not balanced by the 
excitation rate (due to the stochastic forces) in the system. 
Here we attempt to correct this imbalance by "fine-tuning" 
the dissipation rate by an auxiliary thermostat. In order to 
preserve the pairwise conservation of momentum in DPD, this 
auxiliary thermostat is implemented by employing a fluctuat- 
ing dissipation strength, defined by 



7(*) 



2k B T 



-(1 + V (t) At), 



(10) 



where 77 is the thermostat variable. The rate of change of 
X] is proportional to the instantaneous temperature deviation 
f] = C(fcsT — ksT*) where C is a coupling constant, step 
(i) in Table |l[ This first-order differential equation must 
be integrated [step (ii)] simultaneously with the equations of 
motion. In this respect our thermostat resembles the Nose- 
Hoover thermostat familiar from MD simulations. [^9|] Equa- 
tion ( |To| ) can be interpreted as an expansion of the optimal 7 
in terms of At up to the linear order. This ansatz ensures that 
the correct continuum version of DPD is regained for At — > 0. 
Also note that the coupling constant C has to be chosen with 
care. Very small values of C require considerably longer sim- 
ulation times, while too high values may bias the temperature 
distribution as well as the transport coefficients. For the simu- 
lations reported here, we optimized C by studying the charac- 
teristic decay time of (7(^)7(0)). In this manner, we ensured 
that the chosen time scale of the dissipation strength fluctu- 
ations did not interfere with the underlying dynamics of the 
system (in the absence of auxiliary thermostat). 



IV. PERFORMANCE OF INTEGRATORS 
A. Physical quantities studied 

We characterize the integrators by studying a number of 
physical observables. After equilibrating the system, we first 
calculate the average kinetic temperature (ksT), whose con- 
servation is one of the main conditions for reliable simulations 
in the canonical ensemble. Next, we consider the radial dis- 
tribution function g{r), [[30| which is one of the most central 
observables in studies of liquids and solid systems. For the 
ideal gas (model A), the radial distribution function provides 
an excellent test for the integrators, since then g(r) = 1 in the 
continuum limit. Therefore, any deviation from unity has to 
be interpreted as an artifact due to the integration scheme em- 
ployed. For the other models there are no such straightforward 
theoretical predictions. Consequently, we test each model by 
comparing the results of different integrators to one another. 

Artifacts in g(r) are also reflected in the relative isothermal 
compressibility 



1 .ideal 



(11) 
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where K T ieal = (p fc^T*) -1 denotes the compressibility of 
the ideal gas in the continuum limit. For an arbitrary fluid, 
is related to g(r) by 

/>oo 

K T = l + 4np drr 2 [g(r) - 1], (12) 
Jo 

and thus any deviation from Kt = 1 for the ideal gas (model 
A) indicates an integrator-induced artifact. For models B and 
C, kt serves as a measure of integrator-induced artifacts af- 
ter a thorough comparison of results of different integrators 
relative to each other. 

To gauge the underlying problems in the dynamics of the 
system, we consider the tracer diffusion coefficient 

1 N 

D t = }™ g]vt £<[*(*) -*(°)] 2 >' (B) 
i=i 

in which the mean-square displacement ([ri(i) — ^(O)] 2 ) is 
the average squared distance that the tagged particle travels 
during a time interval t. In the long-time limit one obtains 
the tracer diffusion coefficient Dt, which characterizes the 
distance £r> ~ \JDt St travelled by a particle during a long 
time period St. 

Another way to gauge the effects of the numerical inte- 
gration methods on dynamical quantities is to monitor the 
velocity-correlation function 

1 - 

i=l 

which defines the tracer diffusion coefficient through the 
Green-Kubo formula [ (30| 

D T = \ I dt<f>(t). (15) 
Jo 

We note that Eqs. ( pj[ ) and ( p"4| ) are complementary ap- 
proaches for testing the integrators. First, the tracer diffu- 
sion coefficient can easily be measured from simulations via 
Eq. (|l3|), and it provides a way to characterize how possible 
deviations from the true dynamical behavior accumulate to- 
gether. On the other hand, the velocity correlation function 
4>(t) provides relevant information of the short-time dynamics 
of the tagged particle, prior to the region where Eq. ( pjj ) be- 
comes well defined. As an example, the leading term 0(0) 
provides information about temperature conservation, since 
for fluid systems 0(0) = (fcsT)/3m. In addition, since the 
definition of Dt requires 0(t) to decay to zero at long times, 
the decay of <fi(t) can be used to characterize possible short- 
comings in the dynamics of the system. 

B. Results for model A 

First, we discuss the deviations of the observed kinetic tem- 
perature (fcgT) from the DPD-thermostat temperature fcgT*. 
For MD-VV this temperature shift, shown in Fig. [I], is always 



positive and increases monotonically with At. For DPD- 
VV, (kgT) first decreases with increasing At, then exhibits 
a minimum at At w 0.25, and eventually becomes larger 
than IcbT*. The self-consistent approach SC-VV exhibits 
a negative, monotonically increasing temperature shift up to 
At « 0.13, where this scheme becomes unstable at the em- 
ployed particle density. Most importantly, and perhaps most 
surprisingly, we find that the modulus of the temperature devi- 
ation is even larger than for DPD-VV. This finding contrasts 
with the findings of a recent study by Pagonabarraga et al., 
p"2] ] who studied the 2D ideal gas using a self-consistent ver- 
sion of the leap-frog algorithm, and found good temperature 
control for At = 0.06 at p = 0.5. This discrepancy can be 
explained by our observation for the 3D ideal gas that the tem- 
perature shift is in general more pronounced at higher densi- 
ties. A similar effect is found, if the strength of the interac- 
tions is increased. This suggests that temperature deviations 
and other related problems due to large time steps become 
more pronounced when the role of interparticle interactions is 
enhanced. 

In cases where temperature preservation is crucial in cal- 
culating equilibrium quantities, the self-consistent scheme 
with an auxiliary thermostat, SC-Th, is clearly the method of 
choice, as is evident from the results shown in Fig. [l]. For 
this extended-system method, we find that the temperature 
deviations diminish by over two orders of magnitude, with 
a modulus typically of the order of 10 -5 to 10~ 4 . The auxil- 
iary thermostat thus performs very well as long as the iteration 
procedure within the self-consistent scheme remains stable. 

Results for g(r) are shown in Fig. ||. We find that the 
deviations from the ideal gas limit g(r) — 1 are very pro- 
nounced for MD-VV, indicating that even for small time steps 
this integration scheme gives rise to unphysical correlations. 
The performance of DPD-VV is clearly better, while the self- 
consistent scheme SC-VV leads to even smaller deviations. 
For the self-consistent scheme with an auxiliary thermostat 
SC-Th, we found virtually the same results for g{r) as for 
the self-consistent scheme without the thermostat. The results 
for GW and GCC (for a few values of A) integrators were 
approximately the same as those of MD-VV and DPD-VV, 
respectively. For all integrators, the artificial structure in g(r) 
typically becomes more pronounced with increasing time in- 
crement At. It is noteworthy that the bias introduced by the 
self-consistent integrators for At = 0.10 is comparable to that 
introduced by MD-VV already for At = 0.01. 

The relative isothermal compressibilities kt evaluated from 
g(r) are shown in Fig. [3](a). The best performance is found 
for the self-consistent integrators SC-VV and SC-Th, whose 
behavior is essentially similar, and for the DPD-VV whose 
results are almost equally good. In general, the qualita- 
tive behavior of kt reflects our findings for g(r). [ [3l| l The 
magnitude of deviations from kt = 1 is astounding, how- 
ever, and raises serious concern for studies of response func- 
tions such as the compressibility for interacting fluids close 
to phase boundaries. Similarly, the results for tracer diffusion 
in Fig. |3^b) indicate that DPD-VV and the self-consistent ap- 
proach SC-VV work well up to reasonably large time steps, 
while the other integrators were found to perform less well. 



7 



Further studies regarding the decay of the velocity correlation 
function <f>{t) gave similar conclusions, although the size of 
the artifacts in tracer diffusion are best demonstrated by Dt- 
Nevertheless, the decay of velocity correlations in tracer dif- 
fusion is sensitive to the choice of the integrator. 

We now discuss some more general aspects concerning the 
performance of the self-consistent integrator SC-Th. Based 
on our results for (k B T), g(r), and kt, the auxiliary thermo- 
stat performs very well. This provides clearcut evidence that 
the SC-Th scheme is useful for studies of equilibrium quan- 
tities such as the specific heat, which are determined by the 
conservative forces and which are not influenced by the de- 
tails of the dynamics. However, it is less clear whether the 
SC-Th scheme is useful for studies of dynamical quantities. 
To illustrate this point, let us consider the motion of a Brown- 
ian particle as an example. It is characterized by the Langevin 
equation 



M 



dv(t) 
dt 



-Mr)v(t) +F(t), 



(16) 



where M is the mass of the Brownian particle and 77 is the fric- 
tion coefficient which reflects dissipative forces. The remain- 
ing random term F(t) is the driving force due to collisions 
with the solvent particles, whose mass is negligible compared 
to M. In this case, one finds that [B2[| 



D T = 



k B T 
Mr) 



(17) 



which serves to demonstrate that the tracer diffusion of a 
Brownian particle is clearly affected by the dissipation rate. 
Although the motion of DPD particles is not equivalent to 
Brownian motion, the two cases are related. This example 
highlights how any change in the dissipation properties may 
affect diffusion behavior. This problem could arise within the 
SC-Th scheme, since there the strength of dissipation is not 
fixed but it fluctuates in time. Clearly, the significance of this 
issue has to be examined in detail. 

In Fig. ^(a) we show the tracer diffusion coefficient for the 
SC-Th integrator versus the size of the time step At. In the 
inset is shown the average strength of the dissipative force 
(7(t)) as a function of At. The results reveal that Dt con- 
verges to the correct limit at small At, which is expected since 
(7(i)) -> a 2 /2k B T* as At -> 0. For larger time steps, D T 
clearly deviates from the correct behavior. This is due to tem- 
perature deviations ((k B T) < k B T*) within the original SC- 
VV scheme (without an auxiliary thermostat). These devia- 
tions are corrected by the auxiliary thermostat by decreasing 
the average dissipation rate, which in turn increases the diffu- 
sion rate. 

From the discussion above, it is clear that the dissipation 
rate within SC-Th depends on At. Consequently, the trans- 
port properties may not be properly described if temperature 
deviations due to the self-consistent iteration procedure are 
too large. In the present model, this implies that a direct com- 
parison of diffusion properties between SC-Th and other in- 
tegrators is not meaningful. For the purpose of completeness, 
however, let us compare their properties in a slightly modified 



fashion. Instead of comparing the diffusion coefficients them- 
selves, we compare their scaled counterparts Dt^ x / (fc^T). 
This idea is based on an ansatz that tracer diffusion within 
DPD can be written as 



D T 

(k B T) 



(18) 



In Brownian motion, with 7 substituted for 77 in Eq. (|l7|), we 
find the exponent x = 1. For DPD, the situation is differ- 
ent and one finds the behavior of x to be more complex (see 
Sect. |v] and Fig. ^|for further discussion). Detailed studies un- 
der present conditions with the DPD-VV scheme (with small 
At) gave D T /(k B T) ~ I/7 72 . When this dependence on 
the dissipation strength is taken into account, we obtain the 
results shown in Fig. [|(b). Obviously the SC-Th scheme now 
works better, but is nevertheless still less accurate than DPD- 
VV, for example. This finding simply demonstrates that any 
change in dissipation may lead to further changes in the dy- 
namic behavior and should be taken into account in the use 
of auxiliary thermostats. For this reason, we feel that the SC- 
Th scheme is not an ideal approach for studies of dynamical 
quantities by DPD. 

Problems of a similar nature are faced in MD studies with 
the Nose-Hoover thermostat, in which the temperature of the 
system is controlled by a "thermodynamic friction coefficient" 
which is allowed to evolve in time. [|l]| Thus the present prob- 
lem with SC-Th is not speci fic to D PD simulations. Further- 
more, as will be seen in Sect. IV Cj , the SC-Th scheme works 
quite well even for dynamical quantities, when conservative 
interparticle interactions are included in the model. 



C. Results for model B 

We next consider model B, which describes a fluid with 
relatively strong but soft conservative interactions. This situ- 
ation is often met in DPD simulations of polymer dynamics 
and phase separation, among others. Clearly, it is important to 
understand the effects of the integrators on the results in these 
cases. 

As a first and demonstrative topic, we again start by consid- 
ering the deviations of the observed actual temperature (k B T) 
from the desired temperature k B T*. Results shown in Fig. [5] 
for model B reveal that the behavior of the integ rators is very 
similar to that found for model A in Sect. [VB . The largest 
temperature deviations are found for MD-VV and SC-VV, 
and the artifacts due to GCC are almost equally pronounced. 
The performances of DPD-VV and GW are better, while the 
SC-Th scheme is found to be superior to all of them. 

Results for the radial distribution function g(r) resemble 
those for any simple interacting fluid, in this case with a mi- 
nor peak at r w 0.86 r c , and another smaller one around 
r s» 1.55 r c . The radial distribution functions of different 
integrators were essentially similar (data not shown). Thus, 
it is not too surprising that the compressibility data shown 
in Fig. ^(a) do not reveal major differences between differ- 
ent integrators. The results of all integrators are the same to 
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within ±1% for At < 0.01. The differences between dif- 
ferent integrators are very clear only at relatively large time 
steps, where the MD-VV is found to be the poorest and the 
SC-VV the best integration scheme of the ones considered 
here. The results for tracer diffusion in Fig. ||(b) support these 
conclusions. 

The performance of the self-consistent integrator SC-Th 
warrants further attention. We have found that the SC-Th 
provides full temperature conservation for model B. Further- 
more, its results for g(r) and kt are equally good to those 
given by the other integrators, and, finally, even the tracer dif- 
fusion results by the SC-Th are in agreement with results of 
other integration schemes. Thus, for time steps that are not 
too large (say At < 0.01), the self-consistent integrator SC- 
Th seems to provide a promising approach for studies of DPD 
model simulations. These findings contrast with those pre- 



sented in Sect. [V B for the ideal gas. In the present case, the 
differentiating factor is the presence of conservative interac- 
tions. In model B the role of conservative forces is compa- 
rable to the random and dissipative contributions, suggesting 
that the problems in model B due to velocity-dependent dissi- 
pative forces are less pronounced than in model A. Our results 
support this idea. As demonstrated in the inset of Fig. ^(b), 
(7(4)) deviates only slightly (less than 1%) from the desired 
value a 2 /2ksT* at time steps At < 0.01. For larger time in- 
crements, the deviations increase but remain rather small and 
are about 2.5% around At = 0.05. In summary, these re- 
sults demonstrate that, for systems where the role of random 
and dissipative terms is not dominant, the auxiliary thermostat 
not only minimizes temperature deviations, but also provides 
a reasonable approach for calculating dynamical quantities. 



D. Results for model C 

In the preceding models, we have dealt with systems with 
soft interactions. This approach is very suitable for processes 
where the microscopic degrees of freedom do not matter, and 
where one is interested in phenomena at the mesoscopic level. 
However, there are many systems where both microscopic 
and mesoscopic properties are of interest. For example, the 
dynamics of a single polymer chain may be studied using a 
model in which the polymer chain is described in terms of 
Lennard- Jones interactions, while the solvent is treated on a 
mesoscopic level. In this case, the role of DPD would be 
to act as a thermostat and to mediate hydrodynamic interac- 
tions, while the actual interatomic interactions within a poly- 
mer would be described by hard potentials. This approach has 
already proven successful in simulations of a system of small 
amphiphilic molecules, modeled by Lennard-Jones type inter- 
actions in conjunction with the DPD thermostat, [ f33| ] although 
no comparison of the performance of integration schemes was 
reported in that study. 

To clarify the role of integrators in such cases, we exam- 



ine this problem using model C. As described in Sect. [IB 3 



this model uses identical spherical particles, whose pairwise 
conservative interactions are described by a hard repulsive 
Lennard-Jones potential, while the random and dissipative in- 



teractions are soft. Despite its apparent simplicity, this ap- 
proach incorporates the essential aspects required to shed light 
on this issue. 

We focus on two integrators. The MD-VV integrator is 
chosen to represent an approach commonly used in molecular 
dynamics simulations. The performance of MD-VV is then 
compared to that of DPD-VV, which serves as an example of 
integrators designed particularly for DPD. 

We first consider the regime predominated by conservative 
interactions. This is the case for the limit of small a, where 
the role of random and dissipative forces is weak compared to 
that of conservative interactions. The results shown in Fig. ^ 
for the radial distribution function g(r) with a = 1 demon- 
strate that the system is indeed a fluid, and behaves in the ex- 
pected manner. The radial distribution functions for MD-VV 
and DPD-VV are practically indistuingishable, and the same 
holds for the compressibilities extracted from the g(r) data. 
Further studies in this regime revealed that the two integration 
schemes yielded rather similar results for both (ksT) and Dt 
(see results in Fig. ^ as well. Differences between MD-VV 
and DPD-VV are minor at small time steps, but become more 
pronounced as At is increased; around At « 0.01 the de- 
viations are already significant. The temperature conservation 
shows that the artifacts due to MD-VV are stronger than those 
due to DPD-VV. We conclude that in this regime DPD-VV 
performs slightly better than MD-VV. 

The situation becomes more interesting when the random 
and dissipative forces begin to compete with conservative in- 
teractions. The crossover from the regime dominated by con- 
servative Lennard-Jones interactions to the regime dominated 
by dissipative forces takes place around a « 60, as is illus- 
trated in the inset of Fig. ^fb). Right above this threshold, it 
is evident from Fig. ^|(a) that (fc^T) starts to deviate from the 
desired value. In addition, as Fig. ^|(b) reveals, the tracer dif- 
fusion coefficient begins to decrease as soon as a exceeds 60. 
The decrease of Dt simply reflects the fact that the dynamics 
are now governed by random and dissipative forces rather than 
conservative interactions, and thus dissipation slows down the 
motion of DPD particles. Evidently there are similarities with 
Brownian motion, in which Dt ~ I/7. 

The fact that the dynamics in the large-cr regime is con- 
trolled by random and dissipative forces leads us to expect sig- 
nificant quantitative differences between MD-VV and DPD- 
VV, as indeed is observed. First, in this regime the MD-VV 
scheme is less stable than the DPD-VV. Second, tempera- 
ture deviations in the case of DPD-VV are in general less 
pronounced as compared to MD-VV, and the results for the 
tracer diffusion behavior lead to similar conclusions. Thus, 
we conclude that, although the differences between MD-VV 
and DPD-VV are rather small, the DPD-VV method is more 
reliable for simulations in the large-cr regime. 



E. Computational efficiency 

In practice, the choice of an integrator is always a com- 
promise between accuracy and efficiency, which in turn are 
related. Here, we briefly discuss how their mutual outcome 
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can be optimized. 

Based on Tables | and H, it is clear that the efficiency 
of MD-VV, GW, GCC, and DPD-VV is very similar. The 
schemes GCC and DPD-VV require an additional update of 
dissipative forces, but the time it takes is negligible compared 
to the time that is required to update neighbor tables. There- 
fore, these integration schemes are approximately equally ef- 
ficient. The self-consistent approaches, on the other hand, are 
more computer intensive. They are based on an iterative pro- 
cess to find a convergence for dissipative forces and particle 
velocities, an effort which depends on the size of the time step. 
Thus, we focus on a comparison of the efficiency of the self- 
consistent integration schemes to that of DPD-VV. 

Using model A as a test case, we first consider SC-VV. 
As shown in Fig. ^ we find that the SC-VV method requires 
three iterations per integration step to obtain (fcsT) with a 
relative accuracy of 10 -6 at At = 0.01, while 20 iterations 
were necessary at At = 0.10 for the same accuracy. Com- 
pared to DPD-VV, the CPU time per integration step was in- 
creased by a factor of 1.5 for At = 0.10, while it was only 
negligibly higher for the three iterations at At = 0.01. (The 
DPD-VV scheme corresponds to the SC-VV with "zero itera- 
tions".) Fig. ||also shows that the number of iterations needed 
to obtain (fc^T) with a fixed accuracy increases with At, and 
diverges at At « 0.13 where the algorithm becomes unstable 
for the density used here. 

In the case of SC-Th, the CPU time increases by a factor of 
3 to 5 due to the extended simulation times needed to obtain 
the necessary accuracy of (fc^T) with a fluctuating 7. 

These results serve to estimate the computational efficiency 
of the self-consistent integrators in systems where the conser- 
vative force component is very weak. In cases where the role 
of the conservative forces is more pronounced, we expect the 
computational efficiency of the self-consistent schemes to im- 
prove. This is due to the finding that, in models B and C, we 
have seen how temperature deviations in interacting systems 
are smaller than in the ideal gas, and thus a smaller number of 
iteration steps is expected. 

We conclude that the DPD-VV is almost as fast as the MD- 
VV scheme, and the SC-VV scheme is almost as efficient as 
these simple integrators. The SC-Th scheme that includes an 
auxiliary thermostat requires considerably more time. Then it 
is a matter of taste whether the gain in temperature control is 
sufficient to justify the excess in computational effort. 



V. HOW TRACER DIFFUSION RELATES TO THE 
STRENGTH OF THE DISSIPATIVE FORCE 

The tracer diffusion of DPD particles has been the subject 
of various analytical studies. [[7| [34| |35|, ^] Since this topic is 
in part related to the present work (see Sect. IV B ), we wish to 
discuss briefly the relevance of usual approximations made in 
describing the diffusion of DPD particles. 

The descriptions for tracer diffusion of DPD particles are 
usually based on a few reasonable approximations. Most im- 
portantly, the conservative interactions are typically ignored 
and the dynamical correlations between particle displace- 



ments are neglected. Under these circumstances, the system 
is described by the Langevin equation (within the Markovian 
approximation), which yields the tracer diffusion coefficient 



D T 



k B T* 
M7 



(19) 



In practice, this expression describes the diffusion of a Brow- 
nian particle suspended in liquid. In this context, the absence 
of conservative interactions is justified since Brownian motion 
is driven by random forces due to collisions of the Brownian 
particle with the surrounding fluid particles. Neglecting dy- 
namical correlations is also justified, since Brownian motion 
is characterized by a random walk in which case the velocity 
correlation function <p(t) decays exponentially in time, reflect- 
ing the lack of memory effects. 

In models often studied by DPD, the case is rather different, 
however. First, DPD particles move in the presence of simi- 
lar particles, and thus consecutive displacements of the tagged 
particle are likely to be correlated. Second, the conservative 
interactions are not irrelevant. To clarify this issue, i.e., how 
well Eq. ( |l9| ) describes the tracer diffusion of DPD particles, 
we studied the dependence of Dt on the strength of the dis- 
sipative force 7. To this end, we investigated models A and 
B using DPD-VV with a small At (values ranging between 
1 x 10~ 3 - 5 x 10" 3 ). 

The results are presented in Fig. ^. We find that the behav- 
ior of Dt in the two models is very different [Fig. ||(a)]. In 
both cases the power-law dependence Dr/iksT) ~ (l/j) x 
is locally valid, but the exponent x strongly depends on 7 and 
the strength of the conservative force a [see Fig. [j](b)]. 

In the ideal gas {a = 0) the motion of the DPD particles 
is fully governed by the random and dissipative forces, and 
so the exponent x is approximately one at small 7. This be- 
havior is expected, since then the dynamical correlations are 
very weak, as is confirmed by the exponential decay of 4>(t) 
in this regime (data not shown). This is in agreement with 
recent results p5| [36| ] where <p(t) was found to decay expo- 
nentially for small friction. At intermediate values of 7, the 
power-law form of Dt is less clear. The exponent x has a 
minimum around 7=10, and the velocity correlation func- 
tion </>(t) decays algebraically rather than exponentially. Fi- 
nally in the limit of large 7, x tends towards one, which can 
be understood in terms of a large friction force proportional to 
A/7, and therefore this regime mimics the diffusion of DPD 
particles with a large mass. In any case, the decay of <p(t) is 
not exponential, in agreement with analytical predictions by 
Espanol and Serrano. [ |36| l 

In model B with finite conservative interactions, the diffu- 
sion at small 7 is governed by conservative interactions. This 
is best demonstrated in Fig. 0(a), where Dt is only weakly 
dependent on 7 in the limit of small friction. At intermediate 
values of 7, there is a crossover regime in which conservative 
and random forces compete, while at large 7 the diffusion be- 
havior becomes dominated by random and dissipative forces. 
The exponent x varies accordingly, reaching unity only in the 
limit of large 7. 

Our aim in this work is not to focus on the diffusion proper- 
ties of DPD model systems in detail and, thus, we do not con- 
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sider this issue further. Nevertheless, we hope that the present 
results serve to demonstrate that the dynamics in DPD model 
systems are not similar to Brownian motion, and this dissim- 
ilarity is further enhanced by conservative interactions whose 
role can be significant. As regards future studies of transport 
properties of DPD fluid particles, various assumptions made 
in deriving analytical theories should therefore not be taken 
for granted. 



VI. SUMMARY AND DISCUSSION 

Dissipative particle dynamics (DPD) is a very promising 
tool for future large-scale simulations of soft systems. Thus 
far, it has been applied with success to a variety of different 
problems, including studies of pressure profiles inside lipid 
bilayers, [ |To| | phase behavior in surfactant solutions, [[37]] and 
dynamics of polymer chains. [ [3~8| 

Despite its promising nature, DPD has certain practical 
problems that have to be accounted for before extensive use 
in future applications. Many of them are related to the 
coarse-grained nature of the systems studied. As described 
by Espanol and Warren, |Q| the theoretical framework used in 
DPD leads to interparticle interactions that include a dissipa- 
tive term, which depends on the pairwise velocities of DPD 
particles. This implies that for a proper description of the sys- 
tem in time, the dissipative forces and the particle velocities 
should be determined hand in hand in a truly self-consistent 
fashion. As shown in the present work, this issue contains 
various subtle details, but the key point is that the integration 
schemes often used in molecular dynamics simulations cannot 
be used in DPD simulations as such. 

In this work, we have considered this problem through stud- 
ies of three different model systems for a number of inte- 
gration schemes based on the traditional velocity-Verlet ap- 
proach. We have shown that the traditional velocity-Verlet 
scheme gives rise to pronounced artifacts in actual physical 
quantities such as the compressibility and the tracer diffusion 
coefficient. Further studies presented in this work revealed 
that the scale of these artifacts can be greatly reduced by ac- 
counting for the velocity dependence of dissipative forces. 
The simplest approach in this regard is to calculate the dis- 
sipative forces twice during a single time step, and further 
improvements can be obtained if the dissipative forces and 
particle velocities are determined together in a self-consistent 
fashion through a functional iteration process. 

For cases where the remaining temperature deviations need 
to be corrected, we have proposed a self-consistent integrator 
that is coupled to an auxiliary thermostat. We have discussed 
its properties through a detailed analysis in two model sys- 
tems. We have found that this approach works well in the 
case of equilibrium quantities, whose behavior does not de- 
pend on the details of the dynamics. For studies of dynamical 
quantities such as diffusion, however, care must be taken to 
avoid misleading interpretations of the data. Problems may 
appear if (7(i)), extracted from simulations with the auxiliary 
thermostat, deviates significantly from the dissipation strength 
7 determined by the fluctuation-dissipation theorem. In prac- 



tice, this situation is realized if the time step At is relatively 
large and the role of conservative interactions is weak as com- 
pared to random and dissipative forces. However, if care is 
taken and the auxiliary thermostat is used within proper lim- 
its, our results show that it provides an accurate method to 
study DPD models within the NVT ensemble. 

The self-consistent integrator with an auxiliary thermostat 
is an example of a scheme in which the coefficient of the dis- 
sipative force is not constant but fluctuates in time. Conse- 
quently, the average dissipative force strength (j(t)) depends 
on the time increment At. Very recently, den Otter and Clarke 
suggested another approach, Jl4| , [39[ | in which the coefficients 
of the random and dissipative forces depend on the size of 
the time increment At. The results presented in Ref. [14] indi- 
cate that this approach leads to good temperature control [ fto| ] 
compared to the GW scheme, for example, but it remains to be 
shown through thorough tests if this approach is indeed more 
successful in minimizing integrator-induced artifacts than the 
many other schemes suggested previously. 

As noted in the Introduction, DPD can be thought of as 
Brownian dynamics with momentum conservation. Both 
methods are based on coarse graining the underlying micro- 
scopic systems, and in both cases the coarse-grained variables 
are replaced by random noise which is coupled to a dissipa- 
tive friction term. Consequently, one may question whether 
similar integrator-induced problems could be faced in Brow- 
nian dynamics simulations. While we lack direct evidence, 
we feel that the problems in Brownian dynamics (if any) are 
likely less prominent compared to those in DPD. This idea is 
justified by the fact that in Brownian (Langevin) dynamics, 
the velocities of tagged particles are coupled to the dissipa- 
tive force (dili / dt oc —r/Vi) individually for every particle. 
This situation is much easier to deal with compared to that for 
DPD simulations, where the dissipative term includes contri- 
butions from all pairs of particles. A thorough study of this 
topic would be useful. 

In the present work we have examined the performance of 
integration schemes in the well-established description of dis- 
sipative particle dynamics, first suggested by Hoogerbrugge 
and Koelman ||5} and later refined by Espanol and Warren. ^ 
More recently, a number of related schemes have been sug- 
gested to shed more light on the underlying structure of DPD, 
Jl^ ] as well as to generalize the framework of DPD for a num- 
ber of other hydrodynamic cases. pl[ | These approaches are 
numerically more complicated than the DPD considered in 
this work. Studies of the related practical issues would be very 
interesting, although they are beyond the scope of the present 
study. 

We close this work with a brief discussion of the situations, 
in which DPD-specific artifacts due to integration schemes are 
expected. To this end, we first summarize our main findings. 
We have noticed in all three model systems that various inte- 
grators lead to pronounced artifacts in DPD model systems, 
if random and dissipative interactions are strong compared to 
conservative interactions. On the other hand, if the system is 
governed by conservative interactions, then the artifacts have 
been found to be weaker. This suggests that one should use 
conservative interactions that are sufficiently strong to dom- 
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inate the behavior of the model system, and let random and 
dissipative forces act only as a thermostat. Although this ar- 
rangement with dominating conservative forces is feasible in a 
number of cases, [^2[ ^] there are also many systems studied 
by DPD where random and dissipative forces are rather weak 
but comparable to conservative interactions. 

Furthermore there are processes governed by collective ef- 
fects at large particle densities in the high friction limit, which 
based on our work can lead to integrator-induced artifacts. 
Moreover, future aims to examine soft systems of biologi- 
cal molecules in an explicit solvent lead naturally to stud- 
ies of hybrid models, where a microscopic description for 
biomolecules is combined with a coarse-grained description 
for the solvent. In these cases there are both strong conserva- 
tive interactions and relatively weak soft interactions, where 
soft interactions for the solvent are still subject to integrator- 



induced artifacts, and may affect the behavior of the system as 
whole. The overall picture of the role of integration schemes 
in specific model systems is therefore still incomplete, and 
more work is required to resolve these issues, [ pi] Mean- 
while, we emphasize that the artifacts are related to the dis- 
sipative force term, and therefore care should be taken in all 
cases where this term plays an important role. 
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TABLE I: Update scheme for a single integration step (time incre- 
ment At) for various integration schemes in DPD (for acronyms see 
text). For positions and velocities at time t, the updated positions and 
velocities at time t + At are given by the corresponding variables on 
the right-hand side of steps (2) and (4) below, respectively. 
GW(A) : steps (0)-(4), (s) 

MD-VV = GW(A = 1/2) : steps (l)-(4), (s) a 

GCC(A) : steps (0)-(5), (s) 

DPD-VV ee GCC(A = 1/2) : steps (l)-(5), (s) a 

- m + A — (Pf At + P° At + P? \/Al 



m 

m + - — (Pf At + Ff At + VAi 
2 m V 



(0) v° «— 

(1) m «— 

(2) fi < — fi + Vi At 

(3) Calculate Pf{^}, F~P{r 3 ,v°}, ^{r,} 



(4) Vi 



1 1 



2 m 

(5) Calculate F x D {rj, v 3 } 
(s) 6 Calculate k B T~- 



Ff At + F.P At + F t R VAl 



3N-3 



E 



'Sampling step [calculation of temperature kgT , g(r), 
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FIG. 1 : Results for the deviations of the observed temperature (ksT) 
from the desired temperature kaT* = 1 vs. the size of the time step 
At in model A. Results of GW and GCC are for A = 0.65. 



"with substitution of Vj for v° in step (3). 
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TABLE II: Update scheme for the two self-consistent integrators 
without (SC-W) and with (SC-Th) the auxiliary thermostat [steps 
(i)-(iii)]. The self-consistency loop is over steps (4b) and (5) as in- 
dicated. For positions and velocities at time t, the updated positions 
and velocities at time t + At are given by the corresponding vari- 
ables on the right-hand side of steps (2) and (4b) below (after the 
last iteration of the self-consistency loop), respectively. The desired 
temperature is fcsT*. Initialization: rj = 0,7 = a 2 / \2ksT*), and 
UbT is calculated from the initial velocity distribution. 

(i) f) < — C(k B T-k B T*) 

(ii) r] < — T) + 77 At 

(1) v, < — v 4 + i — f Ff At + Ff At + F* VA1 

2 m V 

(2) ft i — fi + Vi At 

(3) Calculate Ff{r 3 }, F t D {r 3 ,v 3 }, F?{r 3 } 

(4a) di < — Vi + i — (jf At + F? VAt ) 
2 m L J 

r (4b) iii < — di + i — At 
2 m 

(5) Calculate ^{fj, Wj} 

JV 

(s) Calculate fc B T= _J?1_ £ i£ , ■•■ 

i — 1 




FIG. 2: Radial distribution functions g(r) vs. At in model A for the 
integration schemes MD-VV, DPD-VV, and SC-VV. Results of GW 
and GCC are almost similar to MD-VV and DPD-VV, respectively, 
and are therefore omitted here. 




FIG. 3: (a) The relative isothermal compressibilities kt vs. At eval- 
uated from g(r) in model A. Ideally one would obtain kt = 1, and 
the deviations from this limit reflect artifacts due to the integration 
procedure, (b) Results for the tracer diffusion coefficient Dt vs. the 
time step At in model A. 




FIG. 4: (a) Results for the tracer diffusion coefficient Dt vs. At 
in model A for the SC-Th integrator, in which 7 is not fixed but 
fluctuates in time. Results of DPD-VV are also given for the purpose 
of comparison. The inset illustrates the dependence of (7(f)) on At 
for the SC-Th integration scheme, (b) The scaled tracer diffusion 
coefficient DtJ x /(feaT) with x = 0.72. 
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FIG. 5: Results for the deviations of the observed temperature (fcsT) 
from the desired temperature fcgT* = 1 vs. the size of the time step 
At in model B. Results of GW and GCC are for A = 0.65. 
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FIG. 6: (a) The relative isothermal compressibilities kt evaluated 
from g(r) in model B. (b) Results for the tracer diffusion coefficient 
D T vs. At in model B. Results of GW and GCC are for A = 0.65. 
The inset illustrates the dependence of {^(t)) on At for the SC- 
Th integration scheme as compared to 7 = 4.5 determined by the 
fluctuation-dissipation theorem. 
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FIG. 7: Results for g(r) in model C with a = 1 using the integrators 
MD-VV and DPD-VV. Results are shown for the density p = 0.1 
with At = 0.01, and for the density p = 0.7 with At = 0.001. The 
results of DPD-VV and MD-VV are essentially identical. 




FIG. 8: Results for (a) the temperature (fcsT) and (b) the tracer 
diffusion coefficient Dt vs. the strength of the random force a in 
model C with p = 0.7. Results are shown for the integration schemes 
MD-VV and DPD-VV with two different time steps. For At = 
0.01 with MD-VV, the system no longer remained stable beyond 
a = 50. To clarify the crossover from the regime dominated by 
conservative interactions to the regime dominated by random and 
dissipative forces, we have shown in the inset of (b) the interactions 
F R (dot-dashed), F D (dashed), and F c (full line) for a = 60. For 
a > 60, the dissipative force is steeper than the conservative one. 
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FIG. 9: Computational efficiency of the self-consistent integration 
scheme SC-VV via functional iteration. Shown here is the aver- 
age number of iterations nu as a function of the employed time step 
At to obtain the desired accuracy. The accuracy is described by the 
modulus of AT /T of the instantaneous temperature, and results are 
shown for AT /T < 10~ 6 (solid diamonds) and AT /T < 10~ 4 
(solid circles). The corresponding CPU time relative to the CPU time 
for plain DPD-VV ("0 iterations") is also shown. 
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FIG. 10: (a) Diffusion results for D T / '(k B T) vs. the strength of the 
dissipative force 7 in DPD simulations for models A (a = 0) and 
B (a = 25). (b) Based on the results in (a) and using the ansatz 
Dt I '{kBT) ~ y~ x , here is shown the resulting exponent 1 as a 
function of 7. The results shown here have been calculated by DPD- 
VV with small At (ranging from 1 x 10 -3 to 5 x 10~ 3 ) such that 
temperature deviations in all cases are very minor. 



